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Abstract 

The mass and velocity distribution of liquid spray has a primary effect on the combustion heat release process. 
This heat release process then affects emissions like Nitrogen Oxides (NO x ) and Carbon Monoxide (CO). 
Computational Fluid Dynamics gives the engineer insight into these processes, but various setup options exist 
(number of droplet groups, initial droplet temperature) for spray initial conditions. This paper studies these spray 
initial condition options using the National Combustion Code (NCC) on a single swirler Lean Direct Injection (LDI) 
flame tube. Using laminar finite rate chemistry, comparisons are made against experimental data for velocity 
measurements, temperature, and emissions (NO x , CO). 


Introduction 

The use of combustion Computational Fluid Dynamics (CFD) in the development of combustion technology has 
been greatly facilitated by the advancements made during the last decade in the areas of combustion modeling, 
numerical simulation, and computing platform. Further development of verification, validation, and uncertainty 
quantification will profoundly impact the reliability and utility of these modeling and simulation tools. Under the 
NASA Fundamental Aeronautics Program, an assessment of existing computational tools for emissions and flow 
field is being carried out. As a first step, the present effort aims at establishing the baseline for prediction methods 
and experimental data for Lean Direct Injection (LDI) (refs. 1, 2, and 3) combustion in confined, swirling flows. 
Combustion codes based on Reynolds Averaged Navier-Stokes (RANS), Very Large Eddy Simulation (VLES), and 
traditional Large Eddy Simulation (LES) will be used; the present paper reports the preliminary investigation using 
the National Combustion Code (NCC). 

The National Combustion Code (NCC) is a state of the art CFD program specifically designed for combustion 
processes. A short summary of the features of the NCC pertaining to this paper are: the use of unstructured grids 
(ref. 4), massively parallel computing — with almost perfectly linear scalability, (refs. 5 and 6) a dynamic wall 
function with the effect of adverse pressure gradient (ref. 7) low Reynolds number wall treatment, (ref. 8), and a 
cubic non-linear k-c turbulence model (refs. 9 and 10), Lagrangian liquid phase spray model (ref. 11), and stiff 
laminar chemistry integration. Recently, viscous low-speed preconditioning (refs. 12 and 13), has been added to 
improve the low-speed convergence of the NCC in viscous regions, and the ability to handle multiple sets of 
periodic boundary conditions has also been added. The combination of these features is usually not available in 
other CFD codes and gives the NCC an advantage when computing recirculating, turbulent, reacting, spray flows. 
Previously, the NCC has undergone extensive validation studies for simple flows (refs. 14 and 15), complex flows 
(ref. 16), NO x emissions prediction performance (ref. 17), and traditional gas turbine combustor/injectors (ref. 18). 

This paper extends the LDI combustion CFD analysis of Davoudzadeh (refs. 19 and 20). Where previous work 
only looked at non-reacting flows, this paper will compare the NCC using a steady- state RANS approach against 
experimental data for a confined, reacting spray flow. This paper will show the sensitivity to the spray initial 
conditions, particularly the number of droplet groups used to define the distribution, and the initial temperature. 
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Figure 2. — Computational grid for the single element LDI combustor. 

I. Flow Modeling 

A. Geometry and Mesh Generation 

The single-element LDI swirler configuration is illustrated in Figure 1 . Each element consists of an air passage 
with an upstream air swirler and a converging-diverging venturi section. The fuel is injected through the center of 
swirler and the fuel tip is at the throat of the venture. The air swirlers have helical, axial vanes with downstream 
vane angles of 60°. There are six vanes with an inside diameter of 9.3 mm and an outside diameter of 22.1 mm. The 
air then dumps into a 50.8 by 50.8 mm combustion section. 

Since the NCC allows for unstructured elements, the grid may be composed of any type and mix of three- 
dimensional elements. However, hexahedral elements were chosen because they are more efficient at filling a 
volume with a smaller number of elements compared to an all tetrahedral grid. Hexahedral elements also allow a 
better calculation of the normal derivatives that are crucial for accurate boundary layer resolution. Approximately 
850,000 elements were used. The “Gridgen” mesh generation software was used to create all the grids used in the 
numerical simulation reported in this paper. 



Figure 1. — Single element LDI injector geometry. 
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B. Air Flow Conditions 


An inlet boundary condition is used, with the air flow speed (20.14 m/s) normal to the inlet face, density 
(1.19 Kg/m 3 ), turbulence intensity level (approximately 5 percent), turbulence mixing length (0.15 m), and static 
temperature (294.28 °K). The inlet pressure in these boundary conditions are calculated, not specified (because they 
are subsonic). The exit boundary condition for single-element specifies static pressure. The operating pressure of the 
combustor is approximately 1 atm, while the equivalence ratio is 0.75. The measured pressure drop (as a percentage 
of P3) during the experiments was measured at 4 percent. 

C. Liquid Spray Conditions 

The spray evaporation process was modeled as a Lagrangian particle, where each particle represents a group of 
actual spray droplets. The droplets are represented as a group rather than an individual droplet because calculating 
each droplet is impractical, as it is too computationally intense. From examining the spray experimental data, a 
Sauter Mean Diameter (SMD, D 32 ) of 32 pm was used for all calculations. The hollow spray cone had a total angle 
of 90.0° and a cone thickness of 14.0°. A mass flow of 4.15xl0” 4 kg/s and an initial velocity of 20.0 m/s was 
prescribed. Simulations were performed with three, nine, and eighteen droplet groups. Figure 3 shows the 
cumulative dropsize distribution (given by Raju (ref. 11)) used as the spray initial condition. Three droplet groups is 
not an adequate representation of the spray distribution; this was used as starting point for calculations. Using nine 
and eighteen droplet groups, the drop size distribution is represented accurately. A well represented (nine droplet 
groups and above) distribution at a 32 pm SMD has a droplet size range from 1 to 95 pm. Ninety-six streams and 32 
rays were used for the spatial distribution for the spray. The Lagrangian solution process uses an “unsteady” spray 
model. Droplet groups are only integrated for a fraction of their lifetime (but restarted at this point for the next 
iteration), rather than completely to steady-state. A time step of 5.0x1 0” 7 was used as the integration time step 
(DTML in NCC terminology), and an integration time of l.OxlO -4 (DGML). New droplet groups were introduced 
every integration time. This process is thought to allow better load balancing on parallel computing platforms, 
possibly allow better coupling with the CFD solver. 


SMD=32, dD=0.2, D max=3 * SM D= 96 micrometers 
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D. Chemistry Model 

Ideally, we would prefer to use detailed chemical kinetic models. There are two problems with this approach: 1) 
Jet-A is a fuel and not a substance, and there are no universally accepted surrogate fuel models for Jet-A; 2) the 
computational costs associated with these models make them impractical when fine computational grids are used. 
Originally, a single-step, global chemistry model was used, shown in Table 1. This model was based on propane 
kinetics (ref. 21), which are close to Jet-A’ s reaction rates. The single-step model allowed an easier start up in the 
solution process, by reducing the computational requirements during the ignition phase. Single-step models do not 
allow emissions calculations, only heat release. Because of this, a reduced ten-step, twelve species model based on 
propane kinetics (refs. 22 and 23) was used, as shown in Table 2. The mechanism was developed by a gradual 
reduction of reaction steps and species using sensitivity techniques. The reduced mechanism also describes the 
formation of Carbon Monoxide and Nitrogen Oxide. However, only one nitrogen-oxide species, NO, has been used 
in the reduced mechanism. NO in the reduced mechanism represents the whole family of nitrogen oxides including 
nitric oxide by Zeldovich (ref. 24) reactions, prompt NO reactions by Fenimore (ref. 25) and nitrogen oxide 
formation through nitrous oxide. 


TABLE L— SINGLE STEP (GLOBAL) CHEMISTRY MODEL 



Reaction 

A, 

mole - cm - sec- K 

n 

E, 

cal/mole 

1 

4 C 12 H 23 + 71 0 2 => 48 C0 2 + 46 H 2 0 
GLO/C 12 H 23 0.10 / 

GLO/0 2 1.65 / 

8.60x10" 

0.00 

3.00xl0 4 


TABLE 2.— REDUCED 10 STEP, 12 SPECIES CHEMISTRY MODEL IN CHEMKIN FORMAT 



Reaction 

A, 

mole - cm - sec - K 

n 

E, 

cal/mole 

1 

4 C 12 H 23 + 47 0 2 => 48 CO + 46 H 2 0 
GLO/C 12 H 23 0.1/ 

GLO/0 2 1.6/ 

1.46xl0 u 

0.00 

3.40xl0 4 

2 

h 2 + o 2 <=> h 2 o + o 

3.98x10" 

1.00 

4.80xl0 4 

3 

H 2 + O <=> H + OH 

3.00xl0 14 

0.00 

6.00xl0 3 

4 

H + 0 2 <=> O + OH 

4.00xl0 14 

0.00 

1.80xl0 4 

5 

CO + OH <=> C0 2 + H 

1.51xlO U7 

1.28 

-7.58xl0 2 

6 

H 2 0 + 0 2 <=> 20 + H 2 0 

3.17xl0 12 

2.00 

1.12xl0 5 

7 

co + h 2 o <=> co 2 + h 2 

5.50xl0 4 

1.28 

-l.OOxlO 3 

8 

N 2 + O <=> N + NO 

l.OOxlO 14 

0.00 

7.50xl0 4 

9 

N + 0 2 <=> NO + O 

631KU) 5 ’ 

1.10 

6.28xl0 3 

10 

N + OH <=> NO + H 

3.80x10" 

0.00 

0.00xl0 u 


E. NCC Computations 

Staging was used in the solution process; cold-flow calculations and initial combustion calculations were 
performed using a single-step chemistry model with Lagrangian spray until a steady state solution was obtained. The 
final stage of CFD calculations was performed by switching from the one- step chemistry model to the reduced 
chemistry model; this was done by changing the input chemistry-parameters of the code. It is important to note that 
no turbulence — chemistry interaction model was used for this case. Based on past experience, this is an adequate 
engineering modeling assumption. 


Results will challenge this assumption. However, the cost of running a complex turbulence-chemistry interaction model is very 
prohibitive, even for supercomputers like NASA’s Columbia. 
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The NCC computations for reacting and non-reacting flow were run in general until the flow residuals were 
reduced three orders of magnitude. The mass flow rates at the boundary conditions were also monitored as a 
convergence criterion. Dissipation (JST type) was set at 0.0 for second order dissipation (c2) and 0.1 for fourth order 
dissipation (c4) (ref. 26). The value of k2, the constant that scales the second order dissipation gradient switch, was 
set at 0.50. Setting the second order dissipation to zero is absolutely necessary to accurately resolving flow features, 
like jets. A CFL number of 1.0 was used. A cubic, non-linear k-c model with a variable C mu coefficient was used. 
This model was selected because of the swirling flow. A dynamic wall function with pressure gradient effects was 
used to model near wall turbulent flow effects. 

Computations were performed on a variety of computer platforms, namely SGI Columbia supercomputer at 
NASA Ames and clusters (Mac OSX and Linux) at NASA Glenn. The Columbia supercomputer was preferred 
because it was considerably faster because of its high speed, low-latency interconnect. This interconnect was 
important because the Lagrangian spray model created a load unbalance in compute nodes. The high speed 
interconnect seemed to mitigate the load imbalance. It takes approximately one week to complete a single LDI 
combustion case using 64 processors. After the baseline run was completed, new runs starting from the baseline case 
took about two days to one week, depending on the number of spray particles used. 

II. Experimental Data 

The experimental data is provided by Jun Cai, S.-M. Jeng, and R. Tacina (ref. 27), and by Yongqiang Fu and 
San-Mou Jeng (ref. 28). Velocity measurements were taken with a two-component Laser Doppler Velocimetry 
(LDV) system, temperature measurements were taken with thermocouples, and emissions data was gathered via an 
isokinetic probe and gas analyzer. Quartz makes up the combustion section. The combustor experiments have 
significant convective and radiative heat losses. The temperature measurements reported are not corrected to 
adiabatic conditions. The heat transfer losses will cause a 200 to 300 °K temperature difference between the 
experimental data and the NCC computations. Experimental droplet measurements are collected with a Phase 
Doppler Particle Analyzer (PDPA). 

III. Results and Discussion 

A. Process Overview 

Table 3 shows the CFD calculations performed to date. In Figure 4, streamlines visualize the flow as it passes 
through the helical swirler and into the combustion chamber. Note the strength of the recirculation zones by the 
track of the streamlines. Contours of axial velocity are also shown. Blue indicates negative axial velocity (and the 
recirculation zones) and red indicates positive axial velocity. The recirculation zone size and strength affects 
combustion stability and emissions. Figure 5 shows the vaporization process. The liquid droplets (small spheres) 
vaporize before combustion occurs. The fuel vapor is shown by the green iso-surface. This fuel vapor mixes further 
with air and then bums. Contours of gas temperature along a centerline cut visualize this combustion process. Blue 
indicates room temperature and red is at the adiabatic flame temperature. Note how the flame shape is bound by the 
fuel vapor region. 


TABLE 3.— SUMMARY OF CFD CALCULATIONS 


Chemistry 

model 

Number of 
droplet groups 

Droplet initial 
temperature, 
°K 

1 Step 

3 

350 

10 Step 

3 

350 

10 Step 

9 

350 

10 Step 

9 

300 

10 Step 

18 

350 
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Figure 4. — Flow visualization using streamlines and contours of axial velocity at an axial mid-plane. 



Figure 5. — Flow visualization of the heat release process using the Lagrangian droplet groups, contours of 
temperature at an axial mid-plane, and iso-surfaces of the fuel mass fraction. 
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B. Flow Comparison 

Figure 6 shows a comparison of axial velocity at the centerline. All of the simulations show general agreement 
with the experimental data, but appear out of phase. The experimental data also shows a dip in axial velocity from 
0 to 10 mm from the combustor face. The calculations do not produce this dip in axial velocity. Figure 7 gives some 
hints to why calculations do not agree. Error bars show the root mean square (RMS) values of axial velocity of the 
experimental data and the calculations. (The RMS values for the CFD calculations are assuming isotropic 
turbulence; RMS values are the same for each velocity and are calculated from the turbulent kinetic energy.) The 
RMS values are very large in the 0 to 10 mm range. Given this high turbulent intensity, the turbulence model (and 
the RANS assumption) will not be able to reproduce this result. Figure 8 shows a spanwise comparison of axial 
velocity. In general, the calculations trend better with the experimental data further downstream. The 3 mm location 
shows a complete disagreement between the experimental data and CFD calculations. This disagreement was also 
show by the centerline comparison. Both nine droplet group simulations match the data better at the 3 mm location, 
but the reason is not understood. Figure 9 shows the comparison with radial velocity for spanwise rakes. The data is 
reasonable matched in all locations. The simulations performed with more droplet groups tend to match the data 
better in the 3 and 9 mm locations. Figure 10 shows a comparison of turbulent kinetic energy along the centerline. 
(The turbulent kinetic energy for the experimental data was created by duplicating the RMS value for the radial 
velocity to get the three values needed to calculate the turbulent kinetic energy. The RMS of axial velocity is the 
other component.) A very interesting result occurs — changing the number of droplet groups drastically changes the 
turbulent kinetic energy produced. The nine droplet group simulations appears to qualitatively match the 
experimental data. Why this exactly occurs is unknown. 1 



Figure 6. — Comparison of axial velocity versus axial distance at the centerline. 


^ While there is no turbulence — chemistry interaction, there is a coupling with the heat release and the turbulence equations. The 
heat release process essentially modifies the expansion and shear the flow is experiencing. This shear then is an input into the 
turbulence equations. We do not claim that a turbulence — chemistry interaction is not needed. 


NAS A/TM— 2008-2 1 5422 


7 



Figure 7. — Comparison of axial velocity with RMS values as error bars 
versus axial distance at the centerline. 
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92 mm Axial Plane 60 mm Axial Plane 28 mm Axial Plane 9 mm Axial Plane 3 mm Axial Plane 



Figure 8. — Comparison of axial velocity versus spanwise distance in the axial Y-Z plane with Z = 0 mm. 
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92mm Axial Plane 60mm Axial Plane 28mm Axial Plane 9mm Axial Plane 3mm Axial Plane 



Distance Across Combustor [mm] 


Figure 9. — Comparison of radial velocity versus spanwise distance in the axial Y-Z plane with Z = 0 mm. 
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Figure 10. — Comparison of turbulent kinetic energy versus axial distance at 
the centerline. 
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Figure 1 1 . — Comparison of temperature versus axial distance at the centerline. 


C. Temperature Comparison 

As a reminder, experimental temperature is not corrected due to heat losses. Comparing temperature along the 
center, Figure 11, as great deal of variation is produced. The one-step model shows a fast reaction rate, while the 
ten-step kinetic model show that process is a function of the number of droplet groups used. While the nine and 
eighteen droplet group distributions do not have a difference in droplet distribution, the eighteen droplet group 
computation shows a much faster reaction rate. In general, the CFD calculations at least qualitatively match the 
shape of the experimental centerline temperature. The more droplet groups used, the better the agreement. Using a 
lower droplet temperature (300 °K) also seems to affect the solution positively. The spanwise rakes, Figure 12 and 
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Figure 13, also compare the CFD calculations against experimental temperature. Generally, the comparison agree 
more downstream. Near the injector, there is a considerable difference. At the 5 mm rake, only the calculations with 
large numbers of droplets agree qualitatively with the experimental data. However, the experimental data shows a 
strong peak temperature at the center of the rake. None of the CFD calculations reproduced this shape accurately. At 
the 10 mm rake, the agreement is better, but only the simulations with adequate droplet distributions (nine and 
eighteen droplet groups), give the qualitative shape. The shape of the peak temperature is still not reproduced 
accurately. 



Distance Across Combustor [mm] 

Figure 12. — Comparison of temperature versus spanwise distance in the axial Y-Z plane 
with Z = 0 mm. 
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Figure 13. — Comparison of temperature versus spanwise distance in the axial Y-Z plane 
with Z = 0 mm. 
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Axial Distance [mm] 

Figure 14. — Comparison of carbon monoxide amounts versus axial 
distance at the centerline. 


D. Emissions Comparison 

Figure 14 compares the amount of CO along the centerline. (Please note: the Y axis is a log scale.) All of the 
10-step simulations reproduce the general amounts of CO. The more droplet groups used, the better the agreement 
with the experimental data. Using a lower initial droplet temperature seems to positively affect the comparison. We 
would like to point out that we have not done any error analysis on this data. Measuring CO is difficult, because it 
requires an isokinetic condition to prevent CO from oxidizing to form C0 2 . Spanwise rakes of CO are show in 
Figure 15 Once again, the more droplet groups used, the better the result. Only the simulations using nine and 
eighteen droplet groups qualitatively agreed with the experimental data, particularly the shapes. While it appears the 
CFD calculations quantitatively agree with the experiments, for the reasons given above, this may not be correct. 
Looking at the NO x amounts along the centerline, Figure 16, we see that the differences increase while proceeding 
downstream. Using more droplet groups does not improve the CFD results. Figure 17 shows the spanwise NO x 
amounts. Close to the injector face, the CFD calculations qualitatively agree with the experimental data. Moving 
downstream from the injector face, the comparison is worse. From this comparison, we believe the ten-step kinetic 
model is not adequate for predicting amounts of NO x . 


NAS A/TM— 2008-2 1 5422 


14 


1 50 mm Axial Plane 80 mm Axial Plane 40 mm Axial Plane 20 mm Axial Plane 



Figure 15. — Comparison of carbon monoxide amounts versus spanwise distance in the 
axial Y-Z plane with Z = 0 mm. 
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NOx 



Figure 16. — Comparison of nitrogen oxides amount versus axial distance at 
the centerline. 
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1 50 mm Axial Plane 80 mm Axial Plane 40 mm Axial Plane 20 mm Axial Plane 



Figure 17. — Comparison of nitrogen oxides amount versus spanwise distance in the axial Y-Z 
plane with Z = 0 mm. 
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Figure 18. — Sauter Mean Diameter (SMD, D32) for an axial distance of 5 mm 
from the combustor face, spanwise rake along the Y = 0 mm. 



Figure 19. — Velocity vectors colored by velocity magnitude in an axial 


slice at the Y=0 mm mid-plane. 

E. Dropsize and Mesh Resolution — Sources of Error 

Looking at Figure 18, we compare computed Sauter Mean Diameter, after the Lagrangian particle is injected into 
the simulation. The CFD and experimental data are measured in essentially the same way — by probes to count each 
droplet size. Details on exactly how this is done will be reported in a later paper. While the overall shape is 
reproduced by the CFD calculations, the dip in experimental droplet size is not captured by the calculations. Also, 
the CFD calculations differ by more than a factor of two. We believe the 32 pm SMD given to the distribution is too 
small. For this case, a 70 pm droplet size would be more appropriate. Figure 19 shows the recirculation zone, via 
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velocity vectors, near the fuel injector tip. The strength of the recirculation zone is large enough to possible cause 
secondary droplet breakup. For a better comparison in the future, secondary droplet breakup should be used. 
Looking at the gaseous fuel mass fraction at the injector tip, Figure 20, there is a considerable amount of evaporation 
that occurs. Computational grid also determines the solution accuracy near the injector tip. Figure 21 shows the 
mesh spacing near the injector tip. Given the amount of change that occurs, the mesh should be refined in this area. 



Figure 20. — Contours of fuel mass fraction (C 12 H 23 ) in an axial slice at 
the Y=0 mm mid-plane. 



Figure 21. — Approximate visualization of the computational mesh in an 
axial slice in at the Y=0 mm mid-plane. 
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IV. Conclusions 


These CFD calculations are a step forward in anchoring CFD codes, like the NCC, to a swirling, reacting spray 
flow. The current results are acceptable, but can and will be improved. The improvements that will be made are: 
using a better reduced Jet-A kinetic mechanism, using a larger Sauter Mean Diameter, adding a secondary droplet 
breakup model, improving the mesh density near the injector tip, and using a consistent turbulence — chemistry 
interaction model. 
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